!=========Calcula período de um pêndulo por integral elíptica==========

PROGRAM exerB2

!======================================================================

IMPLICIT NONE

!========================Declaração de variáveis=======================

      INTEGER :: J
      REAL(8) :: f, h, mid, tetha0, integral, l
      REAL(8), PARAMETER :: g = 9.80665d0

!===========================Leitura dos dados==========================

      WRITE(*,*) "Entre com o comprimento da haste l"
      READ(*,*) l
      WRITE(*,*) "Entre com o ângulo inicial tetha0"
      READ(*,*) tetha0
      WRITE(*,*) "Entre com o h para ser usado na quadratura"
      READ(*,*) h

!==================Faz a integral pela regra de Bode===================

      J = 0
      mid = -tetha0
      integral = 0.0d0
      DO WHILE (mid <= (tetha0-4*h))
            mid = -tetha0 + J*4*h
            !Regra de Bode
            integral = integral + (2*h/45.0d0)*(7*f(mid, tetha0)+32*&
            &f(mid+h, tetha0)+12*f(mid+2*h, tetha0)+32*&
            &f(mid+3*h, tetha0)+7*f(mid+4*h, tetha0))
            J = J + 1
      END DO

!=======================Calcula o período e exibe======================

      WRITE(*,*) "O período do pêndulo é:", DSQRT(2.0d0*l/g)*integral

!======================================================================

END PROGRAM exerB2

!======================================================================

!==============================Funções=================================

REAL(8) FUNCTION f(tetha, tetha0)
      IMPLICIT NONE

      REAL(8) :: tetha, tetha0
      !Define um delta para lidar com as singularidades
      REAL(8), PARAMETER :: delta = 0.001d0
      
      !Faz teste para lidar com possíveis erros
      IF (DCOS(tetha) > DCOS(tetha0)) THEN
            f = 1.0d0/DSQRT(DCOS(tetha)-DCOS(tetha0))
      ELSE IF (DCOS(tetha+delta) > DCOS(tetha0)) THEN
            f = 1.0d0/DSQRT(DCOS(tetha+delta)-DCOS(tetha0))
      ELSE IF (DCOS(tetha-delta) > DCOS(tetha0)) THEN
            f = 1.0d0/DSQRT(DCOS(tetha-delta)-DCOS(tetha0))
      ELSE
            f = 0.0d0
      END IF
END FUNCTION f

!======================================================================
